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Motivation 


Research  Goal: 

Design  a  light  weight  occupant-centric  vehicle  structures  with  high 
levels  of  protection  against  underbody  explosive  events. 

•  A  design  philosophy  that  improve  vehicle  survivability  as  well  as  force 
protection  by  mitigating  soldier  injury  due  to  underbody  IED  &  mine 
blast  is  necessary. 

•  Size  and  weight  would  also  be  paramount  factors  for  Army  vehicles  in 
order  for  faster  transport,  greater  fuel  conservation,  higher  payload 
and  higher  mobility. 

•  Develop  innovative  multilayer  materials  or  structures  to  optimize  the 
dynamic  performance  as  a  mechanism  to  absorb  and  spread  energy 
from  an  impulsive  load. 
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Technical  Approach 


•  Reduced  order  modeling  of  multi-layered  plates  for  rapid  evaluation 
of  the  dynamic  characteristics  for  a  large  numbers  of  alternative 
configurations. 

•  Using  Dynamic  Response  Index  (DRI)  for  injury  as  a  metric  in 
underbody  explosive  events. 

•  Developing  a  meaningful  screening  process  using  the  reduced  order 
models. 

•  Optimizing  the  structural  weight  and  levels  of  protection  of  the 
multilayer  plates  with  a  good  combination  of  materials. 
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Reduced  Order  Model  (ROM) 


Based  on  the  concepts  of  Functional  Graded  Materials  (FGM),  the 
Reverberation  Matrix  Method  (RMM)  is  employed  to  build  a  reduce  order 
model  (ROM)  to  analyze  the  frequency  response  of  an  infinite 
multilayered  medium  under  an  external  excitation. 


The  material  properties  of  each  sublayer,  including  modulus  of  elasticity 
(Et),  Poisson  ratio  (v;),  mass  density  (pt),  loss  factor  (tp),  and  the 
thickness  (hi)  of  sublayer  will  be  assigned. 
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1.  Define  Local  Coordinates  at  Interfaces  and  Boundaries 


/OM) 

x  or  :?:01 
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XN{N- 1) 
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Given  a  solid  medium  with  N  sublayers 
which  has  (N-1 )  interfaces  and  2 
boundaries. 

Define  a  coupled  local  coordinate  systems 
(x,] ,yI])  and  (V(7-1),y/(/-1))  at  each 
interface  -  4  DOFs 

Define  a  regular  local  coordinates  at  top 
boundary  ( x01,y 01)  and  bottom  boundary 

(xw(iv-i)(yw(iv-i))-2DOFs 

Define  p-wave  potential  cp  ( x,y ,  t )  and  s- 
wave  potential  i p(x,y,  t)  in  space  and  time 
domain  for  each  local  coordinate. 
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Reduced  Order  Model  (ROM) 


2.  Double  Fourier  Transform  and  Elastodvnamics  Equations 

•  Wave  potentials  cp  (x,y,  t)  and  ip  {x}y,  t)  satisfy  the  wave  equations: 


v2(P  =  ^  and  V^  =  ^1F 


Where  cp  =  ■  (P-wave  speed)  and  cs  =  :  (SV-wave  speed) 

Apply  double  Fourier  Transform  on  the  wave  equations: 

d2(p  2  ~  n  ,  n?  7  rv 

4-  cczcp  =  0  ;  — |  4-  pzip  =  0 


dy 


dy* 


Where  a  =  I— -  —  k2  :  (p-wave  number  in  y  direction) 


'P 


P  =  -  k2  :  (s-wave  number  in  y  direction) 

And  cp  =  (p(k,y,a)),  \p  =  are  now  the  potentials  in  the 

frequency  (co)  and  wave  number  (k)  domain. 
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Reduced  Order  Model  (ROM) 


2.  Double  Fourier  Transform  and  Elastodvnamics  Equations 


The  solution  of  the  wave  equations  can  be  expressed  as: 

(p(k,y,a))  =  dv{ktcd)e~iay  +  dv(k,od)emy 

\p(k,y,  <x>)  =  ds(k,(jL))e~lPy  +  ds(k,a))el^y 
Where  dp,as  are  the  unknown  arrival  wave  amplitudes; 

dv,ds  are  the  unknown  departure  wave  amplitudes; 
The  displacement  field  u  =  [ux,uy]  and  the 
stress  components  o  =  \oxytoyy\  can  be  calculated  by: 


dcp 
x  dx 


dip  . 
dy  ’ 


chp_  dip. . 

y  dy  dxJ 


•xy 


'yy 


=  ^[2 


d2(p  d2xp  d2xp 


+ 


dxdy  dx2  dy 


-t£]; 


Tr72  ,  n  r d2(P  ,  d2ip 

-AF  <P  +  2/x[^  +  — ]  ; 


in  each  local  coordinates  system. 
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Reduced  Order  Model  (ROM) 


3.  Scattering  Matrices  for  Waves  at  Interfaces  and  Boundaries 


•  The  local  scattering  matrix  SJ  at  interface  J  is  obtained  by  relating 
incident  waves  (arrival)  to  the  transmitted/reflected  waves  (departure). 

•  The  source  wave  vector  sJ  needs  to  be  considered  if  any  wave 
sources  exist  at  interface  J,  where  J  =  1 ,2,..,N. 

•  a  represents  the  “arrival  waves  d  represents  the  “departure  waves  ” 

•  p  represents  the  “P-wave”;  s  represents  the  “SV-wave” 
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Reduced  Order  Model  (ROM) 


3.  Scattering  Matrices  for  Waves  at  Interfaces  and  Boundaries 


+  y 


Interface  J,  y=0 


JU- 1) 


./a+i) 


^  y/Q  +  l) 


The  displacements  and  the  stresses  have 
continuity  relation  at  the  interfaces: 

Ux^-^Qc.a j)  —  uxJV+1\k,a>)  =  0; 
uyyy_1)(fc,  w)  +  uyJ(J+1\k,a))  =  0; 

<*>)  +  axy}(J+1\k,o))  =  0; 

5yyJ<J~1\k,  W)  -  5yyJ(J+1\k,  W)  =  0j 


•  By  substituting  the  elastodynamics  solutions,  the  above  equations  can 

be  expressed  in  a  linear  system  as:  A1  ■  aJ  +  DJ  •  3=  gJ(k,  co)  =  0; 

where  A1  and  D1  are  both  4x4  matrices  which  represent  the  material 
and  physical  properties  of  the  adjacent  layers  of  interface  J. 


•  The  linear  system  can  be  also  written  as :  3  =  SJ  •  a  J  +  sJ] 

where  S1  =  -(Dy)_1-  A1 ,  is  the  local  scattering  matrix  at  interface  J 
sJ  -  -(D^)-1  •  g}  -  0  ,  is  the  local  source  vector  at  interface  J 
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Reduced  Order  Model  (ROM) 


3.  Scattering  Matrices  for  Waves  at  Interfaces  and  Boundaries 


Top  Surface 


Bottom  Surface 


F(a>) 
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4  y(jv+i)jv 


,(JV+l)JV 


-> 


For  the  top  and  bottom  boundaries: 

axy01(k,a) )  =  axyN(N+1fk,a))  -  ayyN(-N+1fk,co)  -  0; 

&yy01(k,a >)  =  F(co); 

5°and  5wnow  degenerate  into  2x2  matrices,  and  the 
source  vector  s°  =£  0  if  loads  apply  on  top  surface. 


•  Construct  a  4Nx4N  global  scattering  matrix  for  entire  multilayered  plate: 
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In  a  compact  notation:  d  -  S  •  a  +  s 
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Reduced  Order  Model  (ROM) 


4.  Phase  Matrix  and  Reverberation  Matrix 


interface  J-l 


x 


0-i)/ 


Layer J 


A 


,/a-D 


interface  J 


The  additional  equation  is  supplemented  by  noting  that 
waves  departing  from  one  side  of  the  layer  become  the 
waves  arriving  at  another  side  of  the  same  layer. 


cLpW  1->  =  eiaJh i  •  dp ^  asJ^  1-)  =  —elPjhi  ■  dj'1  ^  \ 
—  d/(;_1)  =  e~iaJhJ  •  apa~1)J ;  dsJU~1}  =  •  asu~1)J 


S  ° 

Then  all  the  “arrival  waves”  could  be  related  to  the  “departure  waves”  by 
the  global  phase  matrix  P(/i)  as:  a  =  P{h )  •  5* 

where  the  4Nx4N  Global  Phase  Matrix  is: 

0  ...  o  \ 

0  Pi  Oil) 


Local  Phase  Matrix: 

'eiajhj 


Ov)  -  r; 


0 

_eiPjhj 


PQO  = 


0 


P  jv(^v)  0 

0  P  ]\i(hN )/ 
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Reduced  Order  Model  (ROM) 


4.  Phase  Matrix  and  Reverberation  Matrix 

•  The  departure  wave  vectors  d*  and  d  contain  same  elements  but  in 
different  order.  It  can  be  expressed  in  equivalence  through  a  global 
permutation  matrix  U. 

d*  =  U  -d 


where  U  is  a  4Nx4N  block-diagonal  matrix  composed  of  N  same  4x4 
sub-matrix  u\ 


U  = 


0 

u 

0 

0 


0 

0 

u 

0 


,  where  u  = 


0 

0 

0 

1 


•  By  combining  the  global  scattering  matrix  S,  global  phase  matrix  P, 
and  global  permutation  matrix  U,  we  introduce  the  4Nx4N 
reverberation  matrix  R\ 

R(k,(jd)  =  S-P-U 
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Reduced  Order  Model  (ROM) 


5.  Generalized  Ray  Solution  and  Frequency  Response 

•  The  arrival  wave  amplitude:  a  =  P  •  U  •  [I  -  ft]-1*  5 

•  The  departure  wave  amplitude:  d  =  [I  -  s 

Here  the  matrix  [I  -  fl]-1  relates  the  response  of  the  layered  medium 
to  the  excitation  s  in  the  frequency-wave  number  domain. 

•  The  dispersion  relation  for  the  layered  medium  is  found  as: 


det[I  —  R(k}  a))]  =  0; 

Based  on  the  dispersion  relation,  it  is  possible  to  calculate  the 
travelling  wave  number  for  specific  frequencies. 
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Reduced  Order  Model  (ROM) 


5.  Generalized  Rav  Solution  and  Frequency  Response 

•  The  frequency  response  of  displacement  field  could  be  obtained  by 
substituting  a  and  d  into  the  elastodynamics  equation: 

it  *  rr  “I-  Du  •  d 

=  (An-  P  •  d  +  Pu)  •  [/  -  jjF[|T| 

_ V _  v _ 


Transfer  Function 


External  Excitation 
(Input) 


W{k,  y,  co) 


V 


Response  of  the  Plate 
(Output) 


The  matrix  Au  and  Du  are  receiver  matrix 

•  If  we  expand  [/  -  fi]_1  in  a  form  of  Neumann  series: 

[/  -  R]~1=  I  +  R  +  R2  +  ■■■  +  Rm 

Thus  the  response  can  be  evaluated  in  a  summation  form,  where  RM 
means  the  M  times  reflections  and  transmissions  of  the  sources 
waves  in  the  multilayered  plate. 
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Reduced  Order  Model  (ROM) 


r 


5.  Generalized  Rav  Solution  and  Frequency  Response 

•  By  applying  inverse  Fourier  Transform  in  the  wave  number  domain, 
the  frequency  response  could  be  represented  in  space-frequency 
domain: 


Now  the  frequency  response  at  any  locations  in  the  multilayered  plate, 
including  the  top  and  bottom  surfaces  and  all  the  interfaces,  could  be 
calculated. 

•  We  will  validate  the  reduced  order  model  (ROM)  by  comparing  with 
the  spectral  finite  elements  method  (SFEM)  and  numerical  analysis  in 
Nastran  in  the  following  examples. 
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Validation 


1.  Validation  by  Spectral  Finite  Element  Analysis  (SFEA) 


2.  Validation  by  Numerical  Analysis  in  Nastran 

•  Applying  an  impulsive  load 

•  Comparing  the  forced  response  in  frequency  domain 
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Validation  by  SFEA 


1.  Displacement  Field  Approximation 

The  displacement  field  within  multi-layered 
structures  is  approximated  as: 

u(x,y,z,t)  =  N(z)  •  Uelkx~l0)t 

where  u  =  [u  v  w]T  is  the  displacement 
vector,  N(z)  is  the  matrix  of  the  element 
shape  functions,  and  U  is  the  nodal 
displacement  vector. 

For  a  linear  finite  element  n,  shown  to  the  right,  the  displacements  within 
the  element  can  be  expressed  as: 

u(x,y,z,t)  =  N n(z)U n(Jt,  (jS)elkx~l(x)t 

where  the  shape  function  matrix  Nn(z )  =  [(. zn+1  -  z)/hnI  (z  -  zn)/hnI ] 
and  the  amplitude  vector  Un=[Un  Vn  Wn  Un+1  Vn+1  Wn+1]T,  in 
which  I  is  3x3  identity  matrix. 
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2.  SFEA  Formulation 


The  governing  differential  equation  for  free  vibration  of  multi-layered 
structures  can  be  expressed  in  matrix-vector  form  as: 

pit  —  LT  c  •  L  u  =  0 


where  p  is  the  mass  density,  c  is  the  elasticity  constants  matrix,  and  L  is 
the  differential  operator  given  by 


L 


d/dx 

0 

0 


0  0 
0  d/dz 

d/dz  d/dy 


d/dz 

0 

d/dx 


d/dy 

d/dx 

0 


A  standard  one-dimensional  finite  element  formulation  with  the  assumed 
displacement  field  can  be  applied  to  derive  the  following  form  of 
algebraic  eigenvalue  equations: 

[K(k)  -  co2M]U  =  0 
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2.  SFEA  Formulation 


In  the  eigenvalue  equation: 

[K(k)  -  co2M]U  =  0 


the  global  mass  matrix  M  and  the  global  stiffness  matrix  K(/c)  are, 
respectively,  obtained  by  assembling  the  element  mass  matrix  Mn  = 

JQn pnNTnNndz  and  the  element  stiffness  matrix  Kn  =  JQn DTcDdz. 
Here, 


D 


ikNn 

0 

0 


0  0  0  dNJdz  0 

0  0  dNn/dz  0  ikNn 

0  dNn/dz  0  ikNn  0 


Because  the  matrix  D  is  linear  in  k,  it  can  be  shown  that 

K(k)  =  K2k2  +  K±k  +  K0 

Here,  K0,  Klf  and  K2  are  independent  of  the  wavenumber  k. 
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3.  Solution  :  Wavenumber  and  Mode  Shapes 


Based  on  the  eigenvalue  equation,  for  a  given  frequency  co,  is  quadratic 
in  k,  it  can  be  transformed  to  a  linear  eigenvalue  problem  of  the  form: 


[A  -  kB]d  =  0 


where  A  = 


It  is  possible  to  seek  eigenvalues  k  (i.e.  wave  numbers)  and  its 
corresponding  eigenvectors  (i.e.  nodal  displacement  amplitude  vectors) 
U  at  each  frequency  co. 

The  eigenvector  Um  gives  the  deformation  shape  (a.k.a.  mode  shape) 
for  the  /77-th  wave  number  km  at  each  frequency  co.  Eigenvalues  with 
zero  imaginary  part  and  positive  real  part  are  wave  numbers  for 
propagating  waves  and  those  with  non-zero  imaginary  part  are  wave 
numbers  of  evanescent  waves. 
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Validation  by  SFEA 


•  Calculate  and  compare  the  travelling  wave  numbers. 

•  All  tests  are  under  the  frequency  f  =  1 000Hz. 

Case  1:  All  three  layers  are  isotropic  (E  =  200e9  Pa,v  =  0.3, p  =  7800  k9/mi) 


Material 

RMM  Code 

SFEA  Code 

Baseline  (regular  steel) 

20.323  20.323 

Double  Modulus  of  Elasticity 

17.071  17.071 

Half  Modulus  of  Elasticity 

24.204  24.204 

Case  2:  Change  material  properties  in  mid-layer  (top  &  bottom  unchanged) 


Material 

RMM  Code 

SFEA  Code 

Double  E /  Half  E 

19.991/20.525  19.991/20.525 

Double  p/Half  p 

22.118/19.215  22.118/19.215 

v=0.2/v=0.4 

20.337/20.299  20.337/20.299 

Case  3:  Change  material  properties  in  top  &  bottom  (mid-layers  unchanged) 


Material 

RMM  Code 

SFEA  Code 

Double  E /  Half  E 

17.232/23.799  17.232/23.799 

Double  p/Half  p 

22.882/18.576  22.882/18.576 

v=0.2/v=0.4 

20.577/19.947  20.577/19.947 
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Validation  by  Nastran 


A  1m*1m  layered  plate  is  created  in  Nastran  as  another  reference. 


> 


Changing  the 
modulus  of  elasticity 
of  the  middle  layer 


•  The  thickness  for  each  layer  are:  h±  =  h3  =  0.03m,  h2  =  0.04m. 

•  An  impulse  load  is  applied  at  the  center  of  the  top  surface. 

•  Clamped  boundary  condition  applies  on  the  edges  of  the  bottom  surface. 
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Validation  by  Nastran 


•  Capture  the  vertical  displacement  at  the  center  of  bottom  surface. 

•  Compare  maximum  vertical  displacement  obtained  from  ROM/Nastran. 

•  All  the  responses  are  normalized  based  on  the  material  properties  of  the 
first  layer  and  the  magnitude  of  the  excitation  using  u*z  -  logC"**/^1) 

•  Nastran  x  RMM 


Modulus  of  Elasticity  of  middle  layer  [Pa] 
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For  underbody  blast  events,  the  lumbar  and  spinal  injury  would 
be  the  most  threats  to  occupant. 


The  Dynamic  Response  Index  (DRI),  which  is  a  dimensionless 
number  proportional  to  the  maximum  spinal  compression,  will  be 
applied  as  standard  injury  metric  in  this  research. 


Normal  Spine 


The  human  lumbar  would 
be  simplified  as  a  1-DOF 
mass-spring-damper  system. 


Back  view  of  spine 
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UPPER 

BODY 

L>' 
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Dynamic  Response  Index  (DRI) 


The  governing  equation  for  the  simplified  lumbar  model  is: 

d2S  dS  „  d2z 

4-  2  •  %  •  (a)  —  +  (o2  •  8  = 


dLz 


dt2  '  dt  dt 2 

2  -  Time-dependent  shock  acceleration  (excitation  from  the  hull) 
^  -  Damping  Ratio  of  human  lumbar  =  0.224 


uj  -  Natural  Frequency  of  human  lumbar  =  52.9  rad/s 


8  -  Lumbar  Compression 

The  DRI  is  based  on  the  maximum  lumbar  compression  (8max): 

!  DRI  =  ^-5max| 

I _ _ _ j 

The  limiting  DRI  is  17.7  with  a  10%  chance  of  serious  injury. 
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Screening  Metric 


Modulus  of  Elasticity 

Baseline 

E  of  layer  1  (Pa) 

220e9 

E  of  layer  2  (Pa) 

220e9 

E  of  layer  3  (Pa) 

220e9 

E  of  layer  4  (Pa) 

220e9 

E  of  layer  5  (Pa) 

220e9 

Eeff 

220e9 

SAE  INTERNATIONAL 
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The  overall  stiffness  Eeff  of  the 
multilayer  plates  keep  constant 
instead  of  intuitively  changing  the 
stiffness  of  one  certain  layer  to 
change  the  dynamic  response  of 
the  plate.  Therefore,  it  yields  the 
constraint: 

n 

^ *  hi0iai 

i= 1 


Case  1 

Case  2 

Case  3 

300e9 

250e9 

250e9 

250e9 

200e9 

220e9 

150e9 

225e9 

205e9 

250e9 

200e9 

220e9 

300e9 

250e9 

250e9 

220e9 

220e9 

220e9 
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Screening  Metric 


O)2 

The  screening  metric  is  obtained  from  the  value  of  DRI  =  —  8max 

9 


ROM 

Screening  Metric 

Ranking 

Computational  Time 

Baseline 

100% 

4 

10.2s 

Case  1 

81.36% 

1 

10.5s 

Case  2 

94.53% 

3 

10.4s 

Case  3 

94.09% 

2 

10.8s 

Nastran 

Screening  Metric 

Ranking 

Computational  Time 

Baseline 

100% 

4 

2,025s 

Case  1 

87.96% 

1 

2,320s 

Case  2 

97.18% 

3 

2,298s 

Case  3 

95.33% 

2 

2,344s 

•  The  screening  metric  match  well  between  ROM  and  Nastran 

•  The  ranking  order  is  captured  correctly 

•  The  ROM  computes  about  220  times  faster  than  Nastran 
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Design  Optimization 


1.  Optimization  Statement 


An  impulse  load 


I 


An  impulse  load 


N  layers  plate 


STEEL 


I 


mass-spring-damper 

system 

mass-spring-damper 

system 

t0  =  10cm 


Screening  metric 


I 


Screening  metric 


Discrete  Design  Variables:  number  of  layers  (N),  material  selection  (sj) 
Continuous  Design  Variables:  sublayer  thickness  (tt) 

Objective  Function:  min  (Screening  metric ) 

The  screening  metric  is  based  on  the  largest  relative  displacement  experienced  by 
the  spring  mass-damper  system 
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Design  Optimization 


1.  Optimization  Statement 


t0  =  10  cm 


Constraints 

Comments 

lcm  <  tj  <  4.8cm 

Upper  and  lower  bounds  on  the  thickness  of  each  sublayer 

Vw 

9.8 cm  <  2^  ti  <  10.2cm 

Total  thickness  changes  little  to  the  thickness  of  baseline 

2 ‘  Ei*ti>  07  *  ( Esteel *  to) 

Total  surface  stiffness  will  be  kept  in  an  acceptable  range  to  avoid 
reducing  the  plate  rigidity 

2J.  Pi  *  ti  <  (Psteel  *  t0) 

Total  surface  density  should  be  less  than  the  baseline  in  order  to 
avoid  an  increase  in  the  structural  mass 
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Design  Optimization 


1.  Optimization  Statement 

•  Optimization  Solver:  NSGA2  -  type  of  solver  of  the  DS  Toolkit 

•  Generation  size:  10 

•  Population  size  in  each  generation:  200 

•  Material  Selection: 


The  DS  Toolkit  is  used  for  conducting  the  optimization  analysis 
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Design  Optimization 


2.  Optimization  Results 

Convergence  increases  when  the  number  of  generation  increases. 
(The  results  of  the  6  layers  configuration  are  shown  as  an  example) 


1  generation 


2  generation 


3  generation 


4  generation 


5  generation 


6  generation 


7  generation 


8  generation 


9  generation 


10  generation 
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Design  Optimization 


2.  Optimization  Results 

•  3  layers  optimization  (N=3)  •  4  layers  optimization  (N=4) 


steel 

t1=2.32cm 

SiC 

SiC 

t2=4.17cm 

steel 

steel 

steel 

t3=3.60cm 

Ti-Alloy 

t1=2.93cm 

t2=2.64cm 

t3=2.11cm 

t4=2.15cm 


Total  thickness  =  10.1cm 


Total  thickness  =  9.83cm 


Density/unit  area  =  597.31  kg/m  Density/unit  area  =  560.64  kg/m 


•  5  layers  optimization  (N=5) 


t1=1.91cm 

t2=3.20cm 

t3= 1.32cm 
t4=1.78cm 
t5=1.89cm 


Total  thickness  =  10.1cm 


6  layers  optimization  (N=6) 

t1=1 .45cm 
t2=1 .90cm 
t3=1.17cm 
t4=2.29cm 

t5=1.51cm 
t6=1.70cm 

Total  thickness  =  10.0cm 


Density/unit  area  =  520.69  kg/m  Density/unit  area  =  505.27  kg/m 
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3.  Nastran  Validation 


NASTRAN  is  used  to  analyzed  the  optimal  configuration  that  was 
identified  when  a  different  number  of  layers  is  considered.  The  relative 
ranking  of  the  optimal  configurations  that  correspond  to  the  different 
number  of  layers  are  evaluated  and  compared  with  the  ROM. 


In  the  NASTRAN  analysis: 

•  Using  1  m  by  1  m  finite  plate 

•  Applying  the  load  at  the  center  of  the 
lower  surface 

•  Applying  simply-supported  boundary 
condition  at  all  edges  of  the  lower  surface 

•  A  high  structural  damping  boundary  was 
added  to  the  plate  to  reduce  the  reflected 
waves  from  the  edges  of  the  plate  in  the 
NASTRAN  analysis 
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Design  Optimization 


4.  Comparison  Results 


The  screening  metric  is  based  on  the  value  of  the  relative  displacement 
experienced  by  the  single  degree  of  freedom  system: 


120.00% 


100.00% 


80.00% 


60.00% 


40.00% 


20.00% 


0.00% 


ROM  ■  Nastran 
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Q  &  A 


Thank  you! 
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